SALAMANDRA 58(4): 263-274 
30 October 2022 | ISSN 0036-3375 


SALAMANDRA 


German Journal of Herpetology 


Effects of habitat loss on tick load in central populations 
of the Eastern Green Lizard Lacerta viridis 
and its relationship with body condition and population density 


ANA MARÍA PRIETO-RAMiREZ *?, DENNIS RODDER! & KLAUS HENLE? 


53113 Bonn, Germany 


Permoser Str. 15, 04318 Leipzig, Germany 


Manuscript received: 21 September 2021 
Accepted: 27 July 2022 by STEFAN LOTTERS 


Abstract. Habitat loss can increase the susceptibility of individuals to parasitic infections, and hence, parasite load can 
serve as an early warning indicator of stress before the persistence of a population becomes threatened. In this study, we 
tested the effects of patch characteristics, isolation and landscape composition resulting from habitat loss on the tick load 
of individuals from central populations of the Eastern Green Lizard Lacerta viridis. We identified the spatial scale at which 
each landscape composition parameter has the strongest effect and evaluated its effects at this scale. Additionally, we tested 
the relationships between tick load and population density and body condition (BC) to understand possible mechanisms 
that determine tick loads in populations. We found that tick load was not affected by host population density. BC was not 
found to be affected by tick load, but BC did have a negative effect on lizards’ tick loads. The proportion of habitat and 
cropland in the landscape and patch size had positive effects on tick loads, whereas the proportion of urbanized areas, iso- 
lation and perimeter/area ratio had negative effects. We discuss our finding in the context of how the landscape can affect 
tick populations and other potential hosts. We conclude that tick load can be a suitable early warning indicator of nega- 
tive effects of habitat loss, reflecting the susceptibility of lizards to infestation. We suggest that this indicator be included 
in monitoring programs aiming at evaluating the status of populations of L. viridis in modified landscapes, and recom- 
mend that conservation measures be focused on the protection of habitat at broader scales to compensate negative effects 
of cropland and urbanized areas occurring at small scales. 
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of effect, parasites, conservation physiology. 


Introduction 


Habitat loss is the main driver of biodiversity loss, and 
its effects have been mainly studied regarding changes in 
species distribution resulting from the extinction of lo- 
cal populations (ELLis et al. 2012, MAXWELL et al. 2016). 
One of the main aims of conservation biology is to avoid 
the extinction of populations by identifying vulnerable 
populations before they irreversibly decline. Therefore, 
the identification of possible warning indicators of popu- 
lations’ vulnerability is of great importance for improving 
the implementation of monitoring programs in conserva- 
tion. 

Parasitism is a natural ecological regulator of ani- 
mal host populations that has effects on population size 
and demography, reproduction, mating behaviour, mor- 
phology and survival of hosts (POULIN 1999, MOLLER 


2006, BALENGER & ZUK 2014, SANCHEZ et al. 2018). In 
anthropogenically modified landscapes, in which por- 
tions of habitat have been converted into other land cov- 
er types, individuals can be exposed to parasitic infesta- 
tions to a greater extent (Faust et al. 2018). Therefore, in 
these landscapes, the effects of parasitic infestations have 
the potential to be stronger and detrimental to popula- 
tions (Amo et al. 2007, RAHARIVOLOLONA & GANZHORN 
2009). Hence, parasite load has the potential to be useful 
as an early warning indicator of populations’ vulnerabil- 
ity to the negative effects of habitat loss (GILLESPIE et al. 
2005, SANTOS et al. 2018). 

One of the main mechanisms through which an in- 
creased level of parasitism is associated to habitat loss is 
the direct bidirectional link it has with the physiological 
condition of hosts. On the one hand, hosts’ body con- 
dition can be the factor initially being affected by habi- 
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tat loss due to the reduced availability of resources, re- 
duced habitat quality of the remaining habitat (e.g., due 
to edge effects), and increased competition amongst in- 
dividuals (ZANETTE et al. 2000, BUCHER & ENTLING 2011, 
BATTLES et al. 2013). This, in turn, can make individuals 
more vulnerable to parasitic infestations due to immuno- 
suppression (KOSKI & SCOTT 2001, SANCHEZ et al. 2018). 
A higher investment of energy into immune defence is 
required, which consequently reduces mass or fat deposi- 
tions, compromising body condition even more, and in- 
creases susceptibility to parasitic infestations (BOWER et 
al. 2019). On the other hand, an increment in the burden 
of parasites within habitat patches, resulting from para- 
site-spillover effects at newly created ecotones or borders 
after land conversion (FAUST et al. 2018, BORREMANS et 
al. 2019), can occur first; then, this can lead to a reduc- 
tion in the body condition of hosts, which again would 
increase their susceptibility to infestations (SANCHEZ et 
al. 2018). 

The relationship between habitat loss and parasitism 
can also be mediated by host population density, which 
can enhance transmission rates of parasites between in- 
dividuals. This can occur at the initial stages after habi- 
tat loss, when population density transiently increases - a 
phenomenon known as crowding effect (DEBINSKI & HOLT 
2000, TISCHENDORE et al. 2005, KRASNOV et al. 2007) - or 
at later stages, after habitat reduction has occurred, if the 
response of the species to habitat loss still facilitates rel- 
atively high levels of population density (Krasnov et al. 
2002, DEVORE et al. 2020). Even in the case of parasites 
that are not transmitted by direct contact between individ- 
uals, like Ixodes spp. ticks, higher host density can lead to 
increased transmission rates due to elevated indirect con- 
tact rates amongst individuals that share a higher propor- 
tion of space and refuges (OsTFELD et al. 2006, LEU et al. 
2010, FENNER et al. 2011). 

Parasite load and prevalence has been shown to be 
sensitive to different parameters that characterize habitat 
loss. In lizards, it has been demonstrated to increase with 
decreasing habitat quality (Amo et al. 2007, GWIAZDO- 
wicz et al. 2020), and in small and medium-sized mam- 
mals, parasite prevalence and intensity was found to in- 
crease with fragmentation (RAHARIVOLOLONA & GANZ- 
HORN 2009, VAN DER MESCHT 2011). Parameters related 
to landscape configuration and composition have been 
rarely tested, but the few studies in small mammals con- 
ducted to that effect have reported negative correlations 
between habitat availability in the landscape and para- 
site load and prevalence (PUTTKER et al. 2008, THAMM 
et al. 2009, SANTOS et al. 2018). However, although the ef- 
fects of other land cover types on the distribution of para- 
sites have been studied (MAUPIN et al. 1991, TALBOT et 
al. 2019), it has not yet been tested how these land covers 
may affect host parasitism. Moreover, the spatial scale at 
which specific land cover types have the strongest eco- 
logical effects - scale of effect sensu JACKSON & FAHRIG 
(2012) - has been neglected in the research of parasitic 
infestations in transformed landscapes. This is despite 
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the recognized importance of integrating different spatial 
scales to understand the effects of habitat loss (GARMEN- 
DIA et al. 2013, REMM et al. 2017, PRIETO-RAMIREZ et al. 
2020) and of the connection between processes occurring 
at landscape scale and the physiology of individuals (EL- 
LIS et al. 2012). 

In this study, we evaluated the potential usefulness of 
tick loads as an early warning indicator of the vulnerabil- 
ity of populations to habitat loss and fragmentation. To ac- 
complish this, we investigated the effects of patch charac- 
teristic variables, isolation, and landscape composition on 
the tick load (Ixodes ricinus) of individuals of Lacerta viri- 
dis inhabiting a fragmented landscape in central Bulgaria. 
Additionally, in order to gain insights into the possible un- 
derlying mechanisms influencing parasitism, we tested the 
relationship of parasite load with body condition and host 
population density. 

Lacerta viridis is a generalist lizard species with a distri- 
bution range extending from northern Turkey to the cen- 
tral Czech Republic, with the historical core of its distri- 
bution range being located in Bulgaria (MARZAHN et al. 
2016). This species is listed in the European Habitats Direc- 
tive (2007) under Annexure IV, which calls for strict pro- 
tection of the species. Although at a national level, it is the 
most abundant lizard species in Bulgaria (POPGEORGIEV & 
MOLLOV 2005), it does face shrinking habitats due to the 
expansion of land use for agriculture and urbanization in 
the central parts of the country, i.e., on the Thracian pla- 
teau (MOLLOV & GEORGIEV 2015). Earlier studies on liz- 
ards have been focused on the effects of habitat degrada- 
tion and quality on parasite load (Amo et al. 2007, PAFILIS 
et al. 2013, CARBAYO et al. 2019), but the effects of land- 
scape structure on parasite loads in lizards have been stud- 
ied much less (Wu et al. 2019). 

We predicted that parameters that increase resourc- 
es at patch and landscape levels and ensure connectivity 
(patch area, habitat quality, and habitat availability in the 
landscape) have negative effects on the tick load on indi- 
viduals of L. viridis. Conversely, parameters that reduce 
resource availability and connectivity (isolation, patch 
edge, and the proportion of cropland and urbanized ar- 
eas in the landscape) were expected to increase tick load 
(FISCHER & LINDENMAYER 2007, PARDINI et al. 2017, HAT- 
FIELD et al. 2018). Additionally, we expected the amount 
of available habitat in the landscape, patch area, and habi- 
tat quality to reduce the effects of isolation and cropland 
and urbanized areas in the landscape on tick load, and 
edge effects to be intensified by non-habitat land covers 
surrounding patches (HATFIELD et al. 2018). This work 
will contribute to the knowledge on the effects of habitat 
loss on the physiology of L. viridis, the interactions be- 
tween ectoparasites and reptile host species, and the po- 
tential of applying a landscape perspective to investigate 
the physiological status of reptiles in modified landscapes. 
Moreover, it will support evidence-based decisions about 
the usefulness of tick load as an early warning indicator to 
be included in monitoring programs aiming at the protec- 
tion of L. viridis. 


Effects of habitat loss on tick load in Lacerta viridis 


Materials and methods 
Study area 


Our study was conducted in the surroundings of Plovdiv, 
Bulgaria, on the Thracian plateau. This is an alluvial plain 
dominated by the banks of the Maritza River and vari- 
ous tributaries. In this region, which forms part of the 
current and historical core of the distribution range of L. 
viridis (MARZAHN et al. 2016), this lizard inhabits a broad 
range of natural and semi-natural habitat types with high 
vegetation structure, including river banks, shrublands 
and mesophilic mixed forest (MoLLOV 2011). Due to 
urban and agricultural expansion this species faces se- 
vere habitat loss in this region (KAMBOUROVA-IVANOVA 
et al. 2012, MOLLOV & GEORGIEV 2015). Most of the for- 
mer habitat is today configured in fragments of variable 
size that are separated from one another by a matrix of 
uninhabitable land covers. Based on other studies in the 
same area, the process of conversion and fragmentation 
of habitat started in the early 1930, intensifying during 


the 1990° (MOLLOV & GEORGIEV 2015), and the current 
landscape structure might have changed very little, if at 
all, during the last decade (MOLLOV & VELCHEVA 2010, 
MOLLOV 2011). 


Study design and data collection 


In the context of a broader study on occupancy and densi- 
ty in fragmented landscapes we analyzed 42 habitat patch- 
es (PRIETO-RAMIREZ et al. 2018), and based on observa- 
tions in these patches, six patches were selected to develop 
the present study; an additional patch was selected in a 
later survey (Fig. 1, see Appendix 1 for patch coordinates). 
The main criteria for this selection were: 1. patches with 
an encounter rate during the occupancy survey of > 2 
individuals per hour to ensure capturing individuals for 
the present study would be possible within a reasonable 
time period; 2. selected patches should cover the ranges 
of patch area and isolation in the study region; and 3. se- 


Figure 1. Patches where surveys were carried out in the surroundings of Plovdiv, Bulgaria. Map source: Base map ArcGIS desktop 
(ESRI 2015). 
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lected patches should be broadly distributed throughout 
the landscape in such a manner as to be surrounded by 
different land cover types. All selected patches were sepa- 
rated from each other by agricultural landscape, urban ar- 
eas, and/or highways. 

Fieldwork was performed from the beginning of April 
to mid-May in 2015 and from mid-June to the beginning 
of July in 2017. In this region, males become active after 
the winter pause in mid-March, with females following suit 
at the beginning of April. The reproduction season of this 
species lasts until mid- or late May, and both sexes are ac- 
tive and have a unimodal activity pattern during this time 
(GRIMM et al. 2014). As temperature increases with the ad- 
vance of spring and the beginning of summer, both sex- 
es remain active, but their daily activity pattern shifts to 
bimodal with a pause during midday. Thus, survey times 
were adjusted depending on the season: in the spring of 
2015, surveys were performed from 9:00 though 19:00 h, 
and in the summer of 2017, surveys started very early in the 
morning, at 7:00 h, when temperatures were already suit- 
able for lizard activity (~27°C) and continued until 20:00 h 
with a pause of two to three hours during midday, depend- 
ing on temperature. 

Each population was visited on multiple days until 20 
individuals had been captured. Populations 1, 2, 3, and 4 
were sampled only in 2015, whereas populations 5 and 6 
wee sampled in both years, and population 7 was visited 
only in 2017. For each captured individual, sex, age (adult 
or subadult), and the number of ticks (Ixodes ricinus) it 
hosted were recorded. Individuals were weighed with a 
digital scale to the nearest 0.01 g, and snout-vent length 
(SVL) was measured to the nearest 1 mm. Individuals were 
classified by sex and age based on visual inspection of size, 
colour pattern, and reproductive organs in the field. Then, 
we checked for SVL distribution, and our classification 
in the field matched the information provided by Saco- 
nas et al. (2018) for Greek populations of L. viridis about 
minimum size of sexually mature individuals (adult males 
> 80 mm, adult females > 76 mm). 


Calculation of body condition index (BCI) 


We calculated the scaled mass index (SMI) proposed by 
PEIG & GREEN (2009) as a measurement of BCI, which is 
calculated as follows: 


n: | bsma 
scaled mass index: M; = M; El 


where Mand L, are the body mass and linear measurement 
of size of individual 'i'; L, is any value of L within the range 
of values measured. Firstly, a standardized major axis re- 
gression (SMA) of In-transformed body mass on In-trans- 
formed SVL was performed to identify outliers and then 
refitted without outliers to calculate the scale parameter. 
SMI was calculated for all individuals, using the mean of 
SVL (87 mm)as L, and a scale parameter of 3.285. BCI val- 
ues were log-transformed for further analysis. 
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Estimation of population density 


We applied a distance sampling survey with predetermined 
sampling effort to collect data and estimate the densities of 
L. viridis populations (BUCKLAND et al. 1993, ANDERSEN & 
STEIDL 2020). In the context of a broader study, popula- 
tions 1-6 were surveyed to estimate density. In each patch, 
two 1-hour surveys were carried out, each corresponding 
to 1200 m of total walking distance at a speed of 20 m/min, 
which is slow enough to detect lizards. The 1200 m tran- 
sects were split into several units of variable length. The 
number and lengths of transects depended on the different 
habitat types within each patch and the relative area cov- 
ered by each habitat type. This was estimated using satellite 
imagery (see PRIETO-RAMiREZ et al. 2020). Each transect 
was placed within one habitat type only. A strip of 2.5 m 
wide at each side of the transects was scanned to visual- 
ly search for lizards, and the perpendicular distance to the 
transect line was recorded for each detected lizard. 

We estimated densities with the aid of the software Dis- 
tance (THomas et al. 2010). We applied global detection 
probability estimation with stratified density estimation. 
The best-model framework to estimate detection probabil- 
ity was the Conventional simple-distance sampling (CSD) 
engine with a uniform parametric key function and hermit 
polynomial adjustment term, and with empirical estimated 
variance. We selected the model based on the Akaike Infor- 
mation Criterion (AIC), the coefficient of variation (%CV), 
and the Kolmogorov-Smirnov test (K-S test) of goodness 
of fit (AIC = 1167.2, %CV = 4.02, K-S test: estimate = 0.1199, 
p = 0.092). Values of population density used as predictors 
of tick load are presented in Appendix 2. 


Environmental variables 


To assess the impacts of habitat loss on body condition 
and parasite load we evaluated the effects of variables char- 
acterizing patches, habitat configuration, and landscape 
composition. For populations 1-6, values of all variables 
used herein were extracted from a previous study (PRIE- 
TO-RAMIREZ et al. 2020) in which patch digitalization had 
been performed in ArcMap version 10.3.1 (ESRI 2015) to 
calculate patch characteristics and habitat configuration 
variables, and to classify land cover, in order to calculate 
landscape composition variables. For population 7, digitali- 
zation of the patch was performed in this study in the same 
manner in order to calculate patch characteristics and hab- 
itat configuration, but neither habitat quality nor landscape 
composition variables were available for this patch. Geo- 
metric characteristics of each patch include area, perim- 
eter and perimeter/area ratio. Perimeter/area ratio is one of 
the most common metrics to study edge effects following 
habitat loss and fragmentation (FLETCHER et al. 2007), and 
hypotheses regarding edge effects were therefore tested by 
evaluating the effects of this variable. 

To evaluate the influence of habitat quality, we tested 
the effects of vegetation structure. L. viridis is a mesother- 
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mophilic and mesophilic species, and in the study region 
it is highly dependent on enough vegetation cover that will 
offer shade and elevated humidity levels, as well as on ref- 
uges and open spots for basking (MOLLOV & VELCHEVA 
2015, VACHEVA et al. 2020). Thus, it uses habitats with com- 
paratively higher vegetation structure than those used by 
other lizard species in the region (PETROV 2007, VACHE- 
va et al. 2020). Vegetation structure was calculated as a 
‘foliage height diversity’ index (MACARTHUR & MACAR- 
THUR 1961) based on information gathered at microhabi- 
tat scale in several plots in each patch (PRIETO-RAMIREZ 
et al. 2018). 

Habitat configuration was characterized by patch isola- 
tion, calculated as the edge-to-edge Euclidean distance to 
the nearest patch (isolation). For each patch, the isolation 
values were calculated with respect to the factual nearest 
neighbouring patch. Neighbouring patches are not neces- 
sarily included in the present study, but were part of an oc- 
cupancy study performed earlier in the same region (PRIE- 
TO-RAMIREZ et al. 2020). Landscape composition vari- 
ables included the proportion of crop fields and pastures 
(cropland), the proportion of urbanized areas, and the pro- 
portion of habitat surrounding patches. These predictors 
were calculated at different buffer distances around patch- 
es (scales) selected based on dispersal distances reported 
for L. viridis (SAINT-GIRONS & BRADSHAW 1989, MANGIA- 
COTTI et al. 2013, GRIMM et al. 2014,). Selected scales were: 
50, 150, 250, 500, 750 M, 1, 1.5, 2, 2.5 and 3 km (PRIETO- 
Ramírez et al. 2020). The values of variables correspond- 
ing to each patch are summarized in Appendix 3. 


Statistical analyses 


All statistical analyses were performed with the softweare 
R (R Core Team 2018), using the relevant functions (see 
below). In the second year, only very few juveniles that 
were born in that year were captured, and therefore, for all 
subsequent analyses, we classified individuals into two age 
classes: adults and subadults (including juveniles). We ap- 
plied a Levene’s test of homogeneity of variance to test for 
differences between years. No significant differences were 
found (F = 0.013, p = 0.908), and data was therefore com- 
bined for all subsequent analyses. 

Data from all populations were used to test the effects 
of geometric patch characteristics and isolation, but only 
from populations 1-6 to analyze the effects of population 
density, habitat quality, and landscape composition vari- 
ables. Mixed models with patch identity as random effect 
were applied to analyze the data of all traits. Count data 
of tick load was analyzed with GLMM. Count data is usu- 
ally analyzed with a Poisson distribution, but when para- 
site data is over-dispersed, it is recommended that a nega- 
tive binomial distribution be applied (WILSON et al. 1996). 
Hence, we tested for over-dispersion to define the distri- 
bution to be applied. Our over-dispersion test showed that 
our data is not over-dispersed (dispersion = 0.027, p = 
0.884), and therefore, we performed the analysis with a 


Poisson distribution. We tested the effects of sex and age 
by applying mixed models. To evaluate the relationships 
of tick load with BCI and population density, we first ap- 
plied a Pearson correlation to test correlations between 
tick load and BCI, and tick load and population density. 
Then we tested the effect of tick load on BCI by means of 
LMM. Because the relationship between tick load and BCI 
can be bidirectional, the response of tick load to BCI was 
also tested. We also evaluated the response of tick load to 
host population density. Finally, we tested the additive and 
interactive effects of BCI and population density on tick 
load. 

To test the effects of habitat loss on tick load, we first 
identified the scale of effect of each landscape composi- 
tion variable - cropland, urban and habitat - on tick load 
by fitting univariate GLMMs with each of these variables 
at each scale. The scale with the highest R? value was se- 
lected as the scale of effect (MARTIN & FAHRIG 2012) and 
used in subsequent analyses. We built single- and mul- 
ti-scale models. Single-scale models included patch var- 
iables, isolation and landscape composition variables at 
the respective single scale. For multi-scale models, land- 
scape composition predictors were added at their scales 
of effect. We evaluated separately the individual effects 
of each patch and landscape parameter on tick load and 
the effects of the interaction between each parameter and 
sex and age. To test hypotheses regarding the effects of 
interactions between parameters, we constructed models 
with first-order interactions between patch area, perim- 
eter/area, vegetation structure, isolation, proportion of 
habitat, proportion of cropland, and proportion of urban- 
ized areas. Additionally, the effects of interaction terms 
were evaluated in combination with sex or age (second- 
order interaction terms). To avoid collinearity, we tested 
correlations amongst variables using a Spearman rank 
test. Only models including non-correlated variables (r < 
0.60) were tested. 

A large number of multiple comparisons were per- 
formed to evaluate the effects of habitat loss, which by na- 
ture increases the risk of Type I errors. This implies that 
a family-wise error rate (FWER) correction, like the Bon- 
ferroni alpha adjustment, has to be applied (PEDHAZUR & 
SCHMELKIN 1991). However, there are some disadvantages 
associated with the application of Bonferroni corrections. 
It increases the Type II error rate and makes the signifi- 
cance of a single test highly dependent on the number of 
subsequent tests regardless of the effect of the predictor 
itself (PERNEGER 1998, WHITE et al. 2019). Therefore, to 
avoid high Type I error rates while minimizing the short- 
comings of the application of Bonferroni corrections, we 
followed the common null hypothesis principle (MILLER 
1981, RICE 1989), defining ‘family’ as the subgroup of mod- 
els corresponding to a single prediction (BENDER & LANGE 
2001). Thus, a family of models includes the model testing 
the effect of a single patch or landscape parameter, or of 
an interaction between parameters, and two models testing 
the difference of that effect in sex or age classes. Then, the 
calculated alpha threshold (0.05/3) corresponded to 0.016 
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and was applied to all models testing effects of landscape 
and patch parameters. Also, in order to identify the strong- 
er effects, a Cohen’s f test was applied to evaluate effect 
sizes of significant variables or interaction terms (SELYA et 
al. 2012). Results are focused on those variables (or inter- 
action terms) that had medium to high effect sizes (> 0.1) 
(COHEN 1988). 


Results 


A total of 158 lizards were captured from the seven popula- 
tions investigated, being 114 adults and 44 subadults. Sex- 
ing was possible in 145 of these lizards, of which 66 were 
females and 79 were males. Males had a higher tick load 
compared to females (males: range = 0-35, mean=6.89, fe- 
males: range = 0-23, mean = 3.94, t = 6.178, p = 6.5e"°), and 
adults compared to subadults (adults: range = 0-35, mean 
= 6.1, subadults: range = 0-27, mean = 2.93, t = -13.784, 
p < 2e”°). The total number of individuals excluding popu- 
lation 7 was 141, comprising 110 adults and 31 subadults. 
Sexing was possible in 130 of these individuals, of which 65 
were females and 65 males. 


Relationship of tick load to BCI 
and population density 


BCI values ranged from 1.034 to 1.443 with a mean value 
of 1.215. Correlations of tick load with BCI and population 
density were not significant (BCI: r = -0.4, p = 0.354, popu- 
lation density: r = -0.55, p = 0.257). BCI was not affected 
by tick load (B = -0.0008 + 0.0005,t = -1.445, p = 0.151), 
but BCI had a significant negative effect on tick load (8 = 
-1.969 + 0.558, Z = -3.524, p = 0.0004). Density did not sig- 
nificantly affect tick load (B = 0.001 + 0.002, Z = 0.44, p = 
0.66). We did not find additive (BCI: B = 0.2478 + 0.364, 
Z = 0.68, p = 0.497; population density: B = -3.873 + 3.006, 
Z = -1.288, p = 0.198) or interactive effects (B = -6.083 + 
3.634, Z = -1.674, p =0.094) of BCI and population density 
on tick load. 


Effects of habitat loss on tick load 


The scale of effect of landscape composition variables on 
tick load of L. viridis is shown in Figure 2. The proportions 
of habitat and urbanized areas had small scales of effect, 
150 and 250 m, respectively. In contrast, the strongest effect 
of cropland was found at the largest scale, 3000 m. 

Tick load increased with the proportions of habitat, 
cropland and patch area (Figs 3a-c, % of habitat: 8 = 27.921 
+ 5.709, Z = 4.891, p < 0.001, Cohen's f = 0.208, % of crop- 
land: B = 5.419 + 1.468, Z = 3.692, p = 0.0002, Cohen's f? = 
0.138, patch area: B = 3.337 + 1.002, Z = 3.328, p = 0.0008, 
Cohen's f* = 0.104). In contrast, tick load decreased with 
the proportion of urbanized areas (B = -8.1434 + 1.904, Z = 
4.275, p < 0.001, Cohen's f? = 0.154), isolation, and the pe- 
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rimeter/area ratio (Figs 3d-f). In the case of isolation, tick 
load decreased at different rates for each sex (B = -0.682 
+ 0.12, Z = -5.681, p < 0.001, Cohen’s f? = 0.191), with the 
rate being higher for males than for females. The effect 
of isolation was also significantly different between age 
classes (B = 1.211 + 0.227, Z = 5.313, p < 0.001, Cohen's f? = 
0.194), with tick loads decreasing with increasing isolation 
in adults and slightly increasing in subadults. The decreas- 
ing rate of tick load in response to the perimeter/area ratio 
was significantly different between classes, too (B = 1.758 + 
0.386, Z = 4.553, p < 0.001, Cohen's f = 0.164), being higher 
in adults compared to subadults. 

Tick load also increased with vegetation structure at 
low levels of isolation, but its effect decreased with increas- 
ing isolation (Fig. 4A, B = -89.29 + 24.30, Z = -3.674, p = 
0.0002, Cohen’s f = 0.159). Similarly, the effect of vegeta- 
tion structure decreased with increasing perimeter/area 
ratios (Fig. 4b, B = -80.54 + 27.36, Z = -2.944, p = 0.003, 
Cohen's f = 0.117). The effect of urbanized areas was in- 
fluenced by patch area (B = 36.575 + 11.295, Z = 3.238, p = 
0.001, Cohens f? = 0.143), with the tick load increasing 
with the patch area at any proportion of urbanized areas, 
but at a lower rate at high proportions of this land cover for 
both sexes, with this difference being more pronounced in 
males (Fig. 4c). The effect of the perimeter/area ratio was 
influenced by the proportion of urbanized areas (Fig. 4d), 
with the interaction having a negative effect on tick load 
of both males and females, and the tick load being much 
higher in males at low, but not at high proportions of urban 
areas (B = -40.256 + 11.141, Z = -3.613, p = 0.0003, Cohen's 
P = 0.143). 
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Figure 2. Scale of effect of proportion of habitat (Habitat), pro- 
portion of cropland (Cropland), and proportion of urbanized 
areas (Urban) surrounding habitat patches on the tick load of 
individuals of Lacerta viridis. 
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Figure 3. Effects of patch characteristics, isolation and landscape composition variables on the tick load of individuals of Lacerta viridis. 
Coloured dots in a-f represent mean values. Units of untransformed values: patch area in km? and isolation in m. 


Discussion 


We tested the relationship of tick load in individuals of 
L. viridis with population density and body condition, and 
evaluated the effects of patch and landscape parameters re- 
sulting from habitat loss on tick load. In general, we found 
that the latter seems to be related to how ticks experience 


the landscape. We also found a unidirectional relationship 
between tick load and body condition. 

Contrary to our expectation, population density did 
not affect tick load. This might be related to low density 
levels of the studied population. In lizards, tick infesta- 
tion has been found to be related to the number of ticks 
neighbouring individuals host (FENNER et al. 2011). There- 
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fore, with lower population density, the possibility of hav- 
ing close neighbours that influence tick loads will decrease. 
The highest population density in our study (1953.5 ind./ 
km’) was already much lower than the density reported 
for other, less disturbed populations of L. viridis in Bul- 
garia (4535-6000 ind./km?, PopGEORGIEV & MOoLLov 
2005), and similar to that of populations of the closely 
related green lizard species Lacerta bilineata in modified 
landscapes (935-1467 ind./km?) in which habitable patches 
are surrounded by urbanized and agricultural land covers 
(Maura et al. 2011), which suggests that density has prob- 
ably already been affected by habitat loss. 

Also contrary to our prediction, we did neither find a 
significant correlation between tick load and body con- 
dition nor a negative effect of tick load on BCI. However, 
we found significant evidence of a negative effect of body 
condition on the tick load of lizards. Thus, our results sug- 
gest that those lizards having already a low body condition 
index are the ones that are more vulnerable to tick infes- 
tations, but not that tick load compromises body condi- 
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tion. A lower immune response has been linked to lower 
body condition and lower resistance to tick infestations 
in other lacertid lizards, like Podarcis muralis and Lacerta 
agilis (OLSSON et al. 2005, Amo et al. 2006). Tick load has 
been found to have a time-lapsed effect on body condition 
across years in other lizard species (GODFREY et al. 2010), 
which is an aspect we could not test in our study. 

In our study system, the relationship between body con- 
dition and tick load was mediated by characteristics of the 
modified landscape. Body condition was negatively affect- 
ed by the proportion of cropland surrounding patches (un- 
published data), a parameter that we found to have a posi- 
tive effect on tick load. Hence, the plausible underpinning 
mechanism defining the relationship might be a reduction 
of lizards’ physiological status due to negative effects of 
cropland, which in turn renders them more vulnerable to 
tick infections. Moreover, a high proportion of cropland 
around patches can also cause an increment of tick bur- 
den in the landscape due to the presence of other host spe- 
cies. For instance, cropland is known to favour the pres- 
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Figure 4. Effects of interaction terms on the tick loads in individuals of Lacerta viridis. The values of vegetation structure (veg_str), % 
urbanized areas in c), and perimeter/area ratio in d), are the first quantile, mean and third quantile of the data. Units of untransformed 


values: patch area in km? and isolation in m. 
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ence and abundance of rodents (FRAs et al. 2014, BENEDEK 
& SÎRBU 2018, MENDOZA et al. 2019). Rodents are known 
to be ‘bridge hosts’ - host species that cross ecosystem 
boundaries (RENWICK & LAMBIN 2013, BORREMANS et al. 
2019) and are amongst the most likely reservoirs transfer- 
ring Ixodes ricinus ticks from ecosystems such as pastures 
and hedgerows around crops to woodlands (Boyarp et al. 
2008), hence influencing the parasite load in other species 
(BOLANOS-Garcfa et al. 2018). Additionally, some of the 
forested patches visited in this study hold populations of 
deer, a host of Ixodes ricinus ticks whose presence is a key 
determinant of tick burdens in the habitat (GILBERT et al. 
2012, HOFMEESTER et al. 2017). 

Regarding other parameters of habitat loss, we found 
unexpected positive effects of the proportion of habi- 
tat and patch area on tick load. These predictors can also 
increase the tick burden in the landscape and have been 
found to have positive effects on the occupancy and abun- 
dance of ticks (ALLAN et al. 2003, LAWRENCE et al. 2018, 
LEDGER et al. 2019, SIMPSON et al. 2019). Predictors nega- 
tively affecting the tick load of lizards are also related to 
negative effects on tick abundance and lizard population 
density. Contrary to our predictions, tick load was reduced 
by the proportion of urban areas, isolation and perimeter/ 
area ratio. Ticks may be less prevalent in urban wildlife 
due to a lack of habitat that would be suitable for them 
(Murray et al. 2019, DEVORE et al. 2020) and might be 
affected by a lack of habitat connectivity (PUTTKER et al. 
2008). The negative effects of perimeter/area ratios on tick 
loads may correspond to edge effects causing higher turn- 
over and extinction rates in tick populations (SIMPSON et 
al. 2019). 

Literature about the effects of land conversion and habi- 
tat fragmentation on tick loads in individual animals has 
shown that in the absence of a positive relationship be- 
tween host density and parasite intensity, tick load in- 
creases with increasing isolation and decreasing patch size, 
habitat quality, and habitat availability at landscape level 
(RAHARIVOLOLONA & GANZHORN 2009, VAN DER MESCHT 
2011, SANTOS et al. 2018). This is an effect that is expected 
based on the lack of resistance of individuals to infesta- 
tions due to the negative effects of habitat loss and frag- 
mentation on their physiological status (KOSKI & SCOTT 
2001, SANCHEZ et al. 2018). However, in the case of lizard 
species, literature to that effect is scarce, and study results 
are contradictory. For instance, contrary to our results, 
PaFILis et al. (2013) found a higher prevalence and inten- 
sity of tick infestations in habitats with lower vegetation 
structure. Similarly, Gw1azDow!Icz et al. (2020) demon- 
strated that tick load was higher in transformed vs. non- 
transformed habitats in populations of Lacerta agilis and 
Zootoca vivipara, but only in the case of females, whereas 
it was the opposite for males. Finally, similar to our results, 
CARABAYO et al. (2018) found a positive effect of habitat 
quality on tick load in both sexes. This reflects the impor- 
tance of accounting for possible positive effects of patch 
and landscape characteristics on both lizards and ticks, 
which makes it less intuitive to expect low tick load levels 


in less isolated and larger patches, and in landscapes with 
lower habitat availability. 

Finally, regarding the scales of effect of land cover types, 
that of proportion of urbanized areas and the propor- 
tion of suitable habitat on tick load was smaller (150 and 
250 m, respectively) than that of the proportion of crop- 
land (3000 m). This might reflect differences on how di- 
rect the effects of these land covers are on tick abundance, 
and in turn, on lizards tick loads. While habitat and ur- 
banized areas may directly affect tick distribution in the 
landscape, cropland effects would be mediated by its influ- 
ence on other host species. Landscape predictors reflect- 
ing complex species interactions have scales of effect that 
reflect the mobility of all the species involved, which is ex- 
pected to be detected at large rather than localized scales 
(MARTIN & FAHRIG 2012). In this case, the high mobility of 
other hosts, like mammals, through cropland areas might 
be possible to be captured only at broader scales. 

Our study focuses on the effects of landscape, patch, 
population and physiological predictors to understand 
tick load in lizard populations inhabiting a fragmented 
landscape. As discussed above, some of the landscape and 
patch predictors that we analyzed (e.g., habitat availabil- 
ity, vegetation structure) can affect tick populations. How- 
ever, other parameters, like microclimatic conditions and 
patch vegetation composition, can also affect tick distribu- 
tion and abundance (WILLIAMS & WARD 2010, DuMAS et 
al. 2022) and, therefore, potentially also the tick load of liz- 
ards. Integrating parameters affecting ecological processes 
and patterns of tick populations with those included in the 
present study may improve the understanding of host/par- 
asite dynamics in fragmented landscapes, and facilitate a 
better identification of host species susceptibility (BOWER 
et al. 2019). 


Conclusions 


The suitability of tick load as an early warning indicator of 
the negative effects of habitat loss on target reptile species 
depends on its correlation with physiological parameters 
like BC, which we encourage be tested in both directions. 
Based on our results in the studied system, tick load is an 
indicator of an already compromised body condition. Ad- 
ditionally, because tick load can have time-lapsed nega- 
tive effects on the body condition of individuals (GODFREY 
et al. 2010), we recommend this indicator be included in 
studies evaluating the status of individuals in populations 
inhabiting modified landscapes across time. Such long- 
term studies are also necessary to unravel the complex eco- 
logical interactions between ticks and host populations in 
fragmented landscapes. Finally, based on the scales of ef- 
fect obtained for each landscape composition parameter, 
we suggest conservation measures for L. viridis in the stud- 
ied system be focused on compensating the negative effects 
of habitat loss and conversion occurring already at a small 
scale by protecting and increasing habitat availability at a 
much broader scale. 
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